/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include "estep.h"
#include "polytomousrasch.c"
#include "normalize.c"

#include <math.h>
#include <string.h>

static const float log1mepsilon = -0.000010000050000333336;/* log (0.99999) */
static const float logepsilon = -11.512925464970228;       /* log (0.00001) */

float
estep (float            logalpha,
       unsigned int     n_labels,
       const float*     tau,
       const float*     loglambda,
       float*           pz,
       const float*     logpriorz,
       int              clamp,
       const Rating*    ratings,
       unsigned int     n_ratings)
{
  float log_n_labels_m1 = (n_labels < 2) ? 0 : fastlog (n_labels - 1);

  if (clamp < 0 || n_labels < 2)
    {
      memcpy (pz, logpriorz, n_labels * sizeof (pz[0]));
    }
  else
    {
      for (unsigned int j = 0; j < n_labels; ++j)
        {
          if (((unsigned int) clamp) == j)
            {
              pz[j] = log1mepsilon;
            }
          else
            {
              pz[j] = logepsilon - log_n_labels_m1;
            }
        }
    }

  for (unsigned int i = 0; i < n_ratings; ++i)
    {
      float alphalambda = fastexp (logalpha + loglambda[ratings[i].worker]);

      for (unsigned int j = 0; j < n_labels; ++j)
        {
          float dpdbeta;
          float ddpddbeta;

#define taufunc(_k)                                             \
  ({                                                            \
     tau[ratings[i].worker * (n_labels - 1) + _k] + _k + 0.5;   \
  })
          pz[j] += 
            log_polytomous_rasch (ratings[i].rating,
                                  j,
                                  alphalambda,
                                  taufunc,
                                  n_labels - 1,
                                  dpdbeta,
                                  ddpddbeta);
#undef taufunc
        }
    }

  return safe_exponentiate_normalize (pz, n_labels);
}
